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«TITLE OT 
«IDENT /2-007/ ; File: OTSPOWDD.MAR Edit: JCW20 


IT 

DEN 
s enerenconcreenceccccecenscrscqonscceonnscqnenenqqosoeqoeseeeoenceoesooeoeces 
COPYRIGHT (c) 1978, 1980, 1982, 1984 BY 


DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASSACHUSETTS. 
ALL RIGHTS RESERVED. 


THIS SOFTWARE IS FURNISHED UNDER A LICENSE AND MAY BE USED AND COPIED 
ONLY IN ACCORDANCE WITH TERMS OF SUCH LICENSE AND WITH THE 
INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER 
COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY 

HER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY 
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THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE 
AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT 


CORPORATION. 


DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS 
SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL. 
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s FACILITY: Language support library - user callable 
; ABSTRACT: 


1) Double base to double power. 
$} Mh base to double power. 
) Double base to floating power. 


Floating overflow can occur 
Undefined exponentiation can occur if: 
1) neettee base 
2) 0 base and power is 0 or negative. 
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VERSION: 2 
HISTORY: 


; AUTHOR: 
Bob Hanek, 3-Mar-83: Version 2 


; MODIFIED BY: 
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-SBTTL HISTORY ; Detailed current edit history 


| 
| 
Edit history for Version 2 of OTS$POWDD 
5-00) Implemented new algorithm. RNH 3-Mar-83 
-002 The code would overflow when underflow was the proper result. This 
» R2, R6 code, which rounded the 


error arouse from the ADDF5 #*x4DC 
#°X4DC0, R6 was performed the value in 


result so that when —_SUBF 

R6 was greater than R2 comeing RO to be positive when it should have 

pgee negative. The bugfix was to change the test for overflow from 

RO to RO (both of these should be of the same sign). JCW 13-MAY-83. 
2-003 Change INDEX table to be local instead of att Also change 

the ss LEB 3 forma "- INDEX(Rx) to INDEXCRx] to avoid Linker 
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errors. LEB 25-May-1 

2-004 Change remaining table references to use [Rx] instead of (Rx). 
Make Al Moe A2_TABLE local symbols instead of globals. 
LEB 26-A - 8 

e 


2-005 rpangs t 
LEB 29-May 

2-006 Added two ROTL #-3,Rx 
from ‘index*2*3' 
INDEX was_not scaled back to yield values of ‘index’ instead of 

indext2*3" because the mathematics of the code used does need the 

value of index*t2*3 in several computations. JCW 7-Jun-1 

2-007 Corrected a bug involving a SYS_F_FLTOVF_F error during a MULD R6, R2. 
Code was added to see if a MTH overflow message or a zero should be 
returned. JCW 19-Jan-1984 


toyer table reference to use [Rx] instead of (Rx). 
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Rx instructions to scale the value of Rx back 

to ‘index’ before A_TABLECRx] is referenced. The 
| 
| 


DOOOOOOOSCOOCSOoOSS 


Soooooocooo 
SooooooocooSo 
20000 C9 CD CD CD CD 


oo 


12 
a ee SeR0MRAFERECTSION w+ OOUOLE eneclSiON w 1g-sep-19gs 01:52:20 YAR/URS mare YOK ge 


«SBTTL DECLARATIONS 


1 ; INCLUDE FILES: 


We 


i 
2 
0 96 > EXTERNAL SYMBOLS: 
00 98 * -DSABL GBL 
4 99 eEXTRN MTHSSSIGNAL ; Math error routine 
100 eEXTRN MTHSK_FLOOVEMAT i; Floating overflow code 
8 101 ~EXTRN = MTHSK_FLOUNDMAT i Floating underflow code 
: § «EXTRN MTHSK_UNDEXP ; Undefined exponentiation code 
000 104 ; 
00 132 ; MACROS: 
000 1 § : 
000 10 SSFDEF ; Define stack frame symbols 
000 108 ; 
44 1% : EQUATED SYMBOLS: 
00000004 0000 111 ° base =4 ; base input formal - by-value 
0000000C 0000 \1¢ exp = 12 3 exponent input formal - by-value 
00000008 0000 11 fexp = 8 3 exponent when base is floating 
000007FC 0000 114 ACMASK = “M< R2, R3, R4, RS, Rb, R7, RB, RY, R10> 
0000 115 3; register saving mask 


it See: 
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TIONS SEP=1984 11: 
117 ; 
118 ; PSECT DECLARATIONS: 
0000 120 -PSECT _OTSSCODE PIC, SHR,QUAD,EXE ,NOWRT 
! 1 : 3; program section for OTS$ code 
1 j 3; CONSTANTS: 
124 ; 
800 196 
09 ! ; The INDEX table gives the appropriate byte offset into the Al and A2 Tables. 
1 ° 
08 0&8 08 08 08 00 00 00 2 130 INDEX: .BYTE “x00, “x00, “x00, “x08, “x08, “x08, “x08, “x 
18 10 10 13 19 19 13 03 00 131 «BYTE “108: “x10, “x10, *x10, “x48 “AIO, “x10, “x 
20 20 20 18 18 18 18 1 01 13 BYTE “X18, “X18, “X18, “X18, *x18, “x20, “x20, *x20 
8 28 28 28 20 20 ° 4 1 ~BYTE “X20, “X20, “X20, “X20, “X28, “X28, “X28, “x28 
0 30 30 30 ; 0 0 134 -BYTE “X28, “X Be “x30, “X30, “X30, “x30, “x30, “x30 
8 38 38 38 8 0 0028 135 -BYTE “xX g. “x50, “X38, “X38, “X38, “X38, “x38, “x38 
40 40 40 40 40 40 40 38 00 : 136 -BYTE “X38, “X40, “X40, “X40, “X40, “X40, “X40, “X40 
48 48 48 48 48 48 48 4 03 137 -BYTE re fe “X48, “X48, “X48, “X48, “X48, “X48, “X48 
50 50 50 33 23 5 48 0040 138 -BYTE “X48, “x50, “x50, “x50, “x50, “x50, “X50, “x50 
58 58 58 5 : : Boe 139 «BYTE aN23 “x50, “x58, “x58, “x58, “x58, “x58, “x58 
60 60 60 60 60 58 050 140 ~BYTE “x58, “X58, “X58, “X60, “X60, “X60, “X60, “X60 
68 4 68 68 60 60 60 60 005 141 ~BYTE “X60, “X60, “X60, “X60, “X68, “X68, “X68, “X68 
70 8 4 68 68 68 68 68 006 16 ~BYTE “X68, “X68, “X68, “X68, “X68, “X68, “X70, “x70 
70 70 70 70 70 70 4 70 B28 14 «BYTE “x70, “X70, “X70, “X70, “x70, “X70, “X70, “x70 
78 78 78 78 78 78 78 78 «(007 144 -BYTE “X78, “X78, “X78, “X78, “X78, “X78, “X78, “X78 
80 80 80 80 80 78 78 78 th 102 -BYTE “X78, “X78, “X78, *x80, “x80, “X80, “x80, “x80 
0080 147 
“tt 148 «ALIGN QUAD 
080 150 
080 151: For k = 0, 1, .... 16, the k-th entry of the Al table is 2*(k/16) rounded 
0080 13 3 to 56 bits and the k-th entry of the A2 table 2°(k/16) - Al_TABLE(k) 
0080 153 ; rounded to 56 bits. 
0080 154; 
S383 122 arte 
00000000 000 rts o80 139 - QUAD 2100 00090000006 080 
487B67CC AAC rt 0 158 -QUAD “%X487B67CCAAC3408 
PROCS asteaee G95S gp AUB EPC aPC 
BBADSi Bp 5 atte A 18) -QUAD “XBB8A95 a3 ty 98 
Aisa) Fenecoce SUMS 16K AMAR CEA NSaRUEasean 
A14 bats ay 4 a 08 164 -QUAD “XA14 4 ay span 
BEG 3F B84 4 BC 165 - QUAD aABES F904F 3408 
C or 38 oA 3 C 198 eQUAD “x0C fe oA 4 Bp 
ee ght C D 16 -QUAD “Xx ope} gue C 
481151F 24804 33 D 198 eQUAD *X8481151F 24804 cf 
sta tia) 4FC40D E 18 eQUAD “Xx DEBCARS 4FC40D 
Be boES AC bG¢ £0 E 1 9 -QUAD “xX tf EC2AC bee EO 
4 9 DD ¢ EA F 17 eQUA “X24 Ha D ¢ 4OEA 
Atay D4OFS F V6 QUAD “*X ore 4257D40F5 
0004100 0 17 -QUAD *X 000000004160 | 
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FUNCTIONAL DESCRIPTION: 
OTSSPOWOR - DOUBLE result = DOUBLE base ** FLOATING exponent 
OTSSPOWRD = DOUBLE result = FLOATING base ** DOUBLE exponent 
OTSSPOWDD = DOUBLE result = DOUBLE base ** DOUBLE exponent 
The DOUBLE result is given by: 
base exponent result 


0.0 
= 0 Undefined Exponentiation 
<0 Undefined Exponentiation 


ooo Oo OOO 
ow 
3 
~< 


Undefined Exponentiation 
> 8 me * LOG2(base)) 
<0 2* (exp * LOG2(base)) 


VVV A 


Floating Overflow and Underflow can occur. 
Undefined Exponentiation can occur if: 
ase is 0 and exponent is 0 or negative 
2) base is negative 


The basic approach to computing x**y as 2“Cy*log2(x)] is the following: 


eb nab db ss a ts 2s ss 2 ts 2 a a ss 2 bs a > 2) tb 2 5 om 


9.090900 69 00 G9 09 Cd Cd C9 Cd C9. 05.G9 OD Cd CD CD Cd CD Gd Cd Cd OD Cd 09 69 Cd Cd Cd C9 Cd CD Cd Cd CD CD C9 CD CO CD CD CD COCO CN CD CD CDOCECDCRODORCDCD OF 


Step 1: Compute log2(x) to sufficient precision to guarantee an 
accurate final result (see below. 

Step 2: Compute ytlog2(x) to at least the accuracy that log2(x) 
was computed. 

Step 3: Evaluate 2° yetege (xis accurate to the precision of the 
datatype in question. 
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; To determine the accuracy to which boge(x) must be computed to, write 
; prhoacin? as I + h, where I is the integer closest to ytlog2(x5, and 
> h = ytlog2(x) - 1 (Note that ih! =< 1/2.) Then 


2*Cytlog2(x)] = 2°¢1 + h) = (2*1)#(2*h). 


; Since the factor 2°1 can be incorporated into the final result by an integer 
g¢9'5 ton to the exponent field, we can assume that the multiplication by. 
“I incurs no error. Thus the total error in the final result is determined 
3 by how accurately 2*h can be computed. If the final result has p fraction 
ts, we would Like h to have at least p good bits. In fact it would be 
; nice if h had a few extra guard bits, say 4. Consequently, we would Like 
h to be accurate to p + 4 bits. 


MOPOPOROPOPOPoNTy 


et e be the number of bits allocated to the exponent field of the data type 

n question. If I requires more that e bits to represent. then overflow or 
underflow will occur. Therefore if the product y*loge(x) has e + p + 4 good | 
; bits, the final result will be accurate. This requires that log2(x) be 
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computed to at least p + e + 4 bits. 
Since log2(x) must be computed to more bits of precision than is available 


in the base data type, either the next level of precision or multi-precision 
arithmetic must be used. We begin by writing 


=e en+1 
lLog2(x) = lLog2(b) + > cl2ne1)#z', 
A=0" 
Where c(1) = 1, and z* = (2/Ln2)C(z-b)/(z+b)]. Hence 
es 2n+1 
log2(x) = lLog2(b) + 2° + > cl2n+1)#z2' 
nzT~ 
= log2(b) + 2" + plz"). 
Note that if p(z') is computed to P bits, and togett) + 2° is computed 
to pte+4 bits and overhangs p(z‘) by e+4 bits, the required accuracy will 


be achieved. Consequently, the essential tricks, are to pick b such that 
the overhang can be achieved and to compute tog2(b) + z' top +e + 4 bits. 
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CALLING SEQUENCE: 


power.wd.v = OTSSPOWDR (base.rd.v, exponent.rf.v) 
power.wd.v = OTSSPOWRD (base.rf.v, exponent.rd.v) 
power.wd.v = OTSSPOWDD (base.rd.v, exponent.rd.v) 


INPUT PARAMETERS: 
Base and exponent parameters are call by value 


IMPLICIT INPUTS: 
none 

OUTPUT PARAMETERS: 
none 

IMPLICIT OUTPUTS: 
none 


FUNCTIONAL VALUE: 
OTSSPOWDR - DOUBLE base ** FLOATING power 
OTSSPOWRD - FLOATING base ** DOUBLE power 
OTSSPOWDD - DOUBLE base ** DOUBLE power 


; SIDE EFFECTS: 


SIGNALS MTHSK_FLOOVEMAT if floating overflow. 

SIGNALS MTHSK_FLOUNDMAT if ATE | underflow. 

SIGNALS MTH$_ONDEXP (82 = ° UNDEFINED EXPONENTIATION') if 
1) Base is 0 and exponent is 0 or negative 
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2) base is negative 
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GTS spOMDD = DOUBLE PRECI i ai eb 91:3 2:30 as Be o V04-00 Page 
-007 OTSSPOWDD = DOUBLE to DOUBLE giving DOUB 6-SEP-1984 11: MTHRTL SRE SOTSPOWDD.MAR: 1 
O7FC 18 44 ENTRY OTSSPOWDR, ACMASK : 
= ¢ 1A rt SUBL ago. § . ; Allocate 5 longwords on the stack 
Oe 9 1D tg MOVL ; Set number of arguments equal to 2 
04 AE at 9 4 MOVD baie AP), base(SP) ; Move base to stack 
OC AE AC 56 44 CVTFD (AP), exp(SP) ; convert FLOATING exponent to DOUBLE 
OOOOO24C"EF 6€ fa : ig cabs s ie ofs$P WDD ; Call OTS$POWDD 
4 
; 18 
O7FC £ 2 eENTRY OTSSPOWRD, ACMASK g 
SE 14 $3 § 4 51 SUBL 690. $F . ; Allocate 5 longwords on the stack 
6— 02 B 6237 26 MOVL ; Set number of arguments equal to 2 
04 AE 04 AC 6 O23A 5 CVTFD bace(AP). ee ote 3 convert FLOATING base to DOUBLE 
OC AE O8 AC 70 8 3F 54 MOVD fexp (AP) (SP) 3; Move exponent to stack 
QOOO024C'EF 6E FA 0244 355 CALLG (SP), OTSSPOUDD : Call OTS$SPOWDD 
= Be iio. 


OTS! 
1-0 


¢ 13 
TSSPOWDD = DOUBLE PRECISION ** DOUBLE PRECISION p 16-SEP-1 AX/VMS Macro V04-00 Pa 10 
ot 807 OTSSPOWBD = DOUBLE to DOUBLE giving DOUB 6-SEP-19 38 0133 5:63 UMTMRTL. See J0TSP OWDD.MAR;1 pt (7) 
4 5 
O7FC ? ° ; -ENTRY OTSS$POWDD, ACMASK : 
4E 362 ; M RO/R1. If x < 0, or x = 0 and y =< 0, return "UNDEFINED 
: ° : EXPONENT TATION error condition, otherwise attempt to compute x*tty 
4 65 ° 
4 66 GET_BASE: 
50 04 aC 70 4 6 mi MOVD base(AP), RO . 3 RO/R1 <== x 
5 68 COMMON 
1B «14 5 8 BGTR DEFINED ; If x > 0 attempt to compute x*ty 
06 19 BLSS UNDEF INED : Branch to error code for x < 0 
Oc AC (73 7 TSTD exp (AP) ; Check sign of y (Note that x = 0) 
01 #15 4 BLEQ UNDEF INED ; Branch to error condition if y =< 0 
74 ; 
4 3 3M oy wae. continues here, this implies that x = 0 and y > 0. Return 
$8 = 
77: . 
7 


04 


: If processing continues here, this implies that an undefined exponentiation 
3 was attempted. Signal error and return 


| 
| 
| 
| 
RET ; Return 
| 
| 


i eed ed eed FF TY SP PP HMMM EM VAN OOOOOOOOoDOWWWWWD WO LSrorommmmmmmnon 


MEW O CONAN EWN OOD NAMES WWM $$ S ODO NAUSEA CO OONOAU EWN OUOONOULSWN—O”0 CH 


0 
0 
is 
025 
5 
5 
5 
025 
025 
025 
025 7 
025 38 
025 Hy 
025 8 
025 38 
gat 
025 86 UNDEFINED: 
50 8000 ef af Bs? 38 rt ge #*x8000, RO : RO/R1 <== Reserved operand 
7E OO°8F 9A 026 8 MOVZBL pige UNDEXP, =-(SP) ; Put error code on stack 
00000000'GF 01 FB 056 3 CALLS , G*ATHSS$SIGNAL :; Convert error number to 32 bit 
026 9 : condition code and signal error. 
026 39 ; NOTE: Second argument is not re- 
026 39 : quired since there is no JSB entry. 
04 ose * RET ; Return 
$58 33 : 
026 397 ; If process tng continues here will attempt to egapute tty ss ra | yelogeta?s: 
3 egin by determining an such that x = * where =< © ds 
026 9 b ini k and f h th f h 
Oser 400 
54 50 + FFFF807 C 656 40 home #*XFFFFBO7F, RO, RS RG 7 Ne f x) 
» r ‘ $ <-- 2°7 lased ex onent of x 
$4 00004080 BF 3 $597 40 SUBL #*X4080, R4 ; oh Gen 3-78 ia teceenemt of 6% 9) 
50 C2 027 $3 SUBL R4, RO 3 RO <-- iy SF aioe field of x) 
406 ; 
1 407 ; We are now coody ts to compute log2(x). This computation is based on the 
0 8 3 following identity 
5 410 ; 2 Re 1 f-a 
8 411 ; lLog2(2*k*f) = k + log2(a) + ----- wom- 2°(2j+1), where 2 = ----- * 
S 2 3 in 1 ae 2j+1 f+a 
4 414 ; We begin by determining a as b*i, where b = 2°(1/16) and i is 
415 ; between 0 ard 16 inclusive. Rescitinslte i is chosen by table look-up 


OTS 


QTSSPOWDD = DOUBLE PRE ots 
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E to DOUBLE giving DOUB 6-SEP-1984 11:28: MTHRTL.SRCJOTSPOWDD.MAR; 1 


~~ 


in such a fashion as to minimize the magnitude of z. Since log2(a) = i/16 
we may write 


log2(x) = k + 1/16 + zep(z*z). 


FVD te te ee ee ee 


SA 50 =sOFFFFFF8O BF CB #OXFFEFF R 
SA. FD72 CF4A 90 INDEXCR1 R10 <== ¢ 
54 SA CET ADDL3 R : R4, =(SP) SP ==> 2°7e(k + 1/16) 
SA SA COFD BFE OC ROTL #3. R10, R10 R qua | 


ble references like t Line below. 


e Linker will cause an error if 


: R10 

7 R10 

3; SP 

3; R10 

3; R10 will be multiplied by 2*3 by 
; ta e 

: Th 

> 3 

3; 6 Cota 


are used instead of (] for these 
ble references. 


We proceed by computing z = (f-a)/(f+a). In order to insure the accuracy of 
the final result, it is necessary to compute z to at least 68 bits. Since no 
back up data type is available, we must compute z in two parts: z = 21 + 22, 


OOCCCGCCCOCCOCOCOOCOOOOOOCoOOo 


wwowowowowowowowowowowowowowowo 


where z1 is the high 24 bits of z and 22 is the low 56 bits of z. Further, 


1 
1 
1 
1 
1 
1 
1 
1 
9 
F 
8 3 
8 3 
: 3 
0298 > to obtain the desired accuracy it is necessary to work with a = al + a2, 
0298 ; where al and a2 are the high and low 56 bits respectively of ‘a’. We begin 
4 3 computing (in single precision) 
oe 44 : 21 = (f = ald/(f + al) 
0 3B 444 : Note that f-al can be computed exactly in 56 bits, but f+al any require 57 
0298 445 ; bits. The 57 bit can be determined by the exclusive or of the low bits of 
0298 446; f and al. 
38 tt: 
58 FDE3 CF4A 7D 0298 449 ova A1_TABLECR10], R8 3 RB/RO <== al 
7E 51 #59 (CD O29E 450 XORL3 RO, RI, -(SP) 3 SP --> XOR of low bits of al and x 
O2A2 451 : (This will used to determine the 
BSAS re : 7 bit of f+al.) 
52 50 58 61 A2 45 ADCDS =R8, RO, R2 3 R2/R3 <-- f + a1 (rounded) 
50 58 $¢ Q02A6 454 SUBD 28, 3 RO/R1 <== f = al (exact) 
58 50 52 47 O2A9 455 DIVF3 R2, RO, RB 3: RB <-- 21 (single) 
59 00 00 03a $38 MOVL #0, R9 3 R8B/R9 <== 21 (double) 
BO 458 : To compute 22 we note 
+} 459 ; 
B 460 ; z= 21 + 22 = (f = al - a2)/(f + al + a2) 
8 BB 461 ; 
: 186 $ ==> 22 = (f - al - a2)/(f + al + a2) = 21 
6 BO 464 ; Now let v = f + al + a2 = v1 + v2, where v1 and v2 are the high 24 and low 
BS 232 ; 56 bits of v respectively. Then 
8 467 ; 22 = C(f = al = zlevl) = (a2 + ziev2)I/v 
8 469 : We begin by computing v1 and f - al - z1*v1 
BO 471° 
54 52 00 BO 472 MOVL R2, R4 ; R4 <== high longword of f + al 


| 
| 
| 
| 
VAL_LOG2: | 
BICLS s*xFFEFEEBO, RO, R10 <== index to INDEX table 
ROVE, INDEXER10], “a1 je2*3 
-3; RI 
| 
| 
| 
| 


TSSPOWDD = DOUBLE PRECISION ** DOUBLE PRE islon 16-SEP-1984 01: AX/VMS Macro V04-00 Page 1 
oT 889 OTSSPOWDD - DOUBLE to DOUBLE piying DO ub aaa 9133 5388 YATURTL. Se SRE SOTSPOWDD.MAR: 1 . (Fy, 
3 0 00 B3 «473 MOVL #0, R3 3 RG/RS <-- yl 
56 4 6 B6 474 SUBD3 R4, R2, R6 3: R6/R7 <== f + a1 = v1 (exact) 

4 8 4 BA 475 RB, RS 3 RG/RS <== z1*v1 (exact) 

0 4 2 8 $78 SUBD R4, RO 3 RO/R1 <== f = a1 = zitvl (exact) 
$f8 : Compute v2 and a2 + altv2 

54 FEG3 CF4A 7D 480 ° MOVO AZ TABLELR10], _R4 : R4/RS <= a2 
6E FFFEFFFF or CA 481 BICL #°REFFEFFFF, (SP) ; Check if w was rounded 
we 4 § BEQL fs ; Branch if not rounded 
56 FEBD CF 6 4 SUBD TWO_M55, R6 ; Correct for_rounding error (exact) 

6 4 6 484 1$: ADDD R4, R6 3 R6/R7 <-- y 

6 8 64 485 MULD RB. R6 3 R6/R7 <== 21*y2 

6 54 60 $38 ADDD R4, R6 3s R6/R7 <== a2 + z1*v2 
Bs ; Compute 22 

50 56 62 490 ° SUBD RS. RO y RO/R1 <-= (f-al-z1*v1)-(a2-z1*Vv2) 

50 52 66 $3) DIVO R2, RO : RO/R1 <== 22 
19 The next step is to compute log2(x) accurate to at least 68 bits. This is 
t3e accomplished as follows, let 
496 w eon y pe cx x) 
497 (2*7)C 1/16 + z2*p(2z*2z)] 
498 *7e(k + 1/16) + (2*%7)*z2( c2ez*2 + + ¢10*2°10) 

1/16) + } ee ris 
1/16) + + 2' 


® 
?,* wee + d10*#2°*10) 


where 2‘ (2*°7*c0)*z and wi and w2 are the high 24 and low * bits of w 
respectively. Note that the choice of ‘a’ ysed in ‘ee. ee uarantees 
that z* overhangs z'*q(z'*z') by at least 13 bits. Hence, if w is computed 
as wi + w2, the necessary 68 bits of accuracy be be obtained. “00=2%6/ n(2). 


We begin by defining 


Sete Ge Ge Ge Ge Ge Ge Ge Ge Ge Ge Ge Ge Ge Ge Ge Ge Ge Ge Ge Ge Ge Ge Ge 


¢ = hig gp 56 bits of (2*7*c0) 
cl = h 28 bits of (2*7*c0) 
c2 = ieee 56 Ri of (2*7*c0) 
then 
z* = (z1 + 22)*(cl + c2) 
= z1*cl + z21*c2 + 22%. 

54 58 FEBS CF 65 MULD3 C1, RB, R4 3 RG/RS <== c1*z1 
58 FEBB CF 64 MULD C2, R8 3 RB/RO <== c2*21 
50 EFAG CF 64 MULD C, RO : RO/R1 <-- ctz5 

50 58 60 ADDD. RB, RO > RO/R1 <== c#z2 + C221 
58 50 54 61 ADDD3 =R4, RO, RB 3 RB/RY <== 2° 


We proceed by letting 
wi = high 24 bits of 2*7#(k + i/16) + z1*c1 
w2" = (C2*7#(k # 1/16) * z1*c1 = wl) + 21*c2) + 22*c. 


TH HAH HHH HMMM MMMMMMMMMMMMMMMmMmMMMmMmmMMMMMmMmMmooOVTVVTIMOOOONO 
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DPVLVIVIVIVIVIVPUPUSVDIVIVSIVSVSUSUSVSVSVSUSVSVSVSVSTSISISISIS 
ROMPINPIMPINIMYNININ 2 2 OO SH MH OOOOOOCOOCCO NO 
WDONAOUNE WIN @ QWDNAULS WN $$ O OONAULS WN O00 


and 


g 
parry 
+ 

=< 
~m 


sso 


4 5 
7 0 
52 6 A 
53 0 
54 2 
6— 50 4 


50 58 
FE93 CF 04 
50 


unum 
mosOoc 


54 OC AC 


52 


53 
53 : 


8 
“4 


3 FT vee Nd 
2 
4 


7? 

8 
2 

5 7F 
6 


=| 
: 


F 
7 
F 
F 
4 
& 
F 
2 
6 
6 


> 


oO 
= 
“oO 
v 
= 
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Sete Se Ge Ge Ge Ge Se Ge 
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; yl and y2 are the hig 
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==> 2*7*(k + 1/16) + 2° = 
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wi + w2'. 


==> w = (C2*7#(k + 1/16) + 2°] + z'*qQ(z'*z') 
w 


= 
= wl + we, 


where w2 = w2" + 2'tq(z'*z') 
CVTLF 4(SP), R10 
ADDF 3 R10, R4, RE 
SUBF3 R16, R6, R2 


MULD3 = RB 


, RB, RO 
POLYD RO, #LOGLEN-1, LOGTAB 


R8, 
ADDD (SP), RO 


+ we’ + 2'8q(z'*2') 


R10 <== 2°7(k + 1/16) 

R6 <== 2*7(k+i/16) + 21*¢1 in single 

sped <== wl 

R2 <== bits of z1*c1 included in wi 
= -C2(k+i/16)+z21%¢1 = wi) 

Convert R2 to double 

C2*7(k+i/16)4z1%c1 = wi] + c1*z1 

(SP) ==> w2' 


RO/R1 <== z'*2' 

RO/R1 <== q(z'*z") 
RO/R1 <== 2'#Q(z'*2z') 
RO/R1 <== w2 


We now calculate yetopeta? = (yl+y2)*(witw2) = ylewl + y2ewl + ytw2, where 


MOVQ exp(AP), R4 


and low 28 bits of y respectively. 


3; R4&/RS <W- y 


Test for the possibility of overflow in the computation of yry!. 
This will occur if the exponent of y plus the exponent of wl is greater 


; than 127 


EXTZV «#7, #8, R4, R2 
SUBW2 #*Xx80, R2 


EXTZV #7, #8, R6, R3 

SUBW2 #*K80, R3 

ADDW : 

CMPU | #*X7F, R3 

BLSS Y_TIMES_W1_OVER 
MULD RS, R 

MOVL. R4, R 

BICL3 #*XFFFFOFFF, RS, R3 
SUBD 2, R4 

MULD = R6. RZ 


/ L " 
MULD R6, R4 


unbiased exp of wl 
unbiased exp of wity ; 
largest unbiased exp possible is 127 


RO/R1 <-- yewe 

<-- high longword of y1 
R2/R3 <== yl 
R4/RS5 <== y2 
R2/R3 <== ylewl 
R4/R5 <== y2ewl 


; The next step in computing 2“*Cy*log2(x)] is to write y*log2(x) as 


—o 
an 


$e 


13 | 
TSSPOWDD = DOUBLE PRECISION ** DOUBLE PRE Fsi0N 16-SEP-1984 01:57: AX/VMS Macro Vv04-00 Pa 14 | 
ts OTSSPOWDD = DOUBLE to DOUBLE oletan DouB ae ets 9 91:38:69 EMTHRTL. SREIOTSPOUDD MAR; 1 - (7) | 
C 73 
: i ; ytlog2(x) = I + j/16 + g/16, 
C : where 1 is an integer, is an integer between 0 and 15 inclusive, and | 
: : : g is a fraction in'the nterval c-172, 172) get 
¢ 93 ° | 
56 52 000040C0 BF 41 C 94 ADDF3 #*X4DC0, R2, R6 : 3*2°S is used in this truncation process | 
re 32 5 Sere 8 gees tehe +: eats a hy 
3 could occur e number is n 
56 000040C0 BF 42 0564 39 SUBF —#X4DCO, RE RG <o= D*2C1 $ 9716) in single ™ 
90 D 68 98 MOVL #0, R $ Rg/R? <== 2°7(] + he: in double 
2 6 $f 3 SUBD R6, R2 3; R2/R3 <== those bits of wityl not 
0371 600 3 included in 2*7(1 + j/16) 
26 52 60 0371 601 ADDD R2, R4 ; 
9 54 60 8 74 one ADOD R4, RO ; RO/R1 <--g 
5 38 49 77 ~=©60 CVTFW R6, R7 3 R7 <== 2*78(1 + j/16) in integer 
5 1D 037A 604 BVS EXCEPTION_1 : Branch if {1} is too Large 
bare 606 ; 
037C 607 ; We can now compute 
037C 608 ; 
ba5¢ 29? ; xtty = 2*Cytlog2(x)] = 2*CI + j/16 + g/16) 
o37¢ 611 : = (2*1)*CAe(B+1)] = 2*1*CA + A®B), where 
037C gi : A = 2°(j/16) is obtained by table look-up and B = 2*(g/16) - 1 is obtained 
037C 614 ; by a Min/Max approximation. 
Osre ele 
FES6 CF eT Be rs Gre sir POLYD RO, #EXPLEN-1, EXPTAB 3; RO <-- B = 2*(g/16) - 1. (Chebyshev 
Baas 618 : polynomial approx. of degree 5) 
52 57 FFFFFF8O 8F CB 038 619 BICL3 #°XFFFFFF8O, R7, R2 é <-=- 2°3 * index into Al_TABLE 
a 52 FD 8F 9C O38A 620 ROTL #-3, R2 3 R2 <== index into A1_TABLE 
0 FCEC CF42 64 O38F 621 MULD  A1_TABLECR2], RO > RO/R1 <== A®B 
28 FD6E CF4 60 395 6 § ADDD A2_TABLECR2J, RO 3; RO/R1 <== AtB + A2 
5 FCEO CF4 60 0398 6 ADDD Al_TABLECR2), RO 3 RO/R1 <-- 2*C(j+g)/16)] = (A*B+A2)+A1 
57 O7F 8F AA O3A1 624 BICW #°R7F, R 3 R7 = 2478] 
0 57 Ad O5A6 625 ADDW = R7,_-RO > RO <== 2*18#2°C(j+g)/16) 
OO7F 8F 50 B11 OQO3A9 626 CMPW RO, #*X7F 3; test for over/underflow 
07 «15 8 ar 627 BLEQ EXCEPTION_2 3 see what exception is if nee or = 0 
04 peat ° 3 RETURN: RET 3 otherwise return result in RO 
1 6303 
351 631 ; Handlers for software detected over/underflow conditions follow 
1 3 : , 
1 634 EXCEPTION_1: 
53 1 635 TSTF R6 : if big ARG > 0 goto overflow 
1D 3618 6 § BGEQ OVER F handler, otherwise go to 
S ee ° , BRB UNDER s underflow handler 
7 639 EXCEPTION 2: 
57 B5 7 640 TSTW R7 ; test sign of I; if 1 <0 
17 «+18 9 oe) BGEQ OVER ; go to overflow handler 
$28 
| 
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low 

® Cw 


BB 644 ; ~ would have caused a hardware gt floating o erf error. If y<0, 
BB 645 ; then we should return a result of 0 since result = 2*(y*(witw2)). Note, 
B 976 3 y can not be zero. 
BB 647 ; 
BB 908 
BB 649 Y_TIMES_W1_OVER: 
54 73 BB 650 TSTD R4 : if y <0 no overflow is needed 
13. (14 BO 651 BGTR OVER : overflow 
i 
: O22 ; Underflow; if user has FU set, signal error. Always return 0.0 
+ 
QO3BF 657 UNDER: 
50 7C O38F 658 CLRQ RO : RO/R1 <-- 0 
0B 04 AD «(06 COE C1 659 BBC #6, SFSW_SAVE_PSW(FP), 2g 
C6 660 ; has user enabled yg A under flow? 
7E OO°8F 9A C6 661 Atty pag FLOUNDMAT, -(SP) ; Put underflow code on st 
00000000 ' GF 01 FB CA 666 CAL es ATHS$SIGNAL ; Convert to MTH$_ FLOUNDMAT “(32-bit 
03D1 66 3; VAX=-11 exception code) and 
0301 664 3 signal condition 
04 03D1 665 2$: RET 3; return 
03D 666 
03D 667 ; 
b30 o38 ; Signal floating overflow, return reserved operand, -0.0 
03D 670 
03D 671 OVER: : else process for overflow 
7E OO°SF QA O03D2 ore MOVZBL #MTHSK_FLOOVEMAT, -(SP) ; Put overflow code on stack 
50 01 OF 79 0306 67 ASHQ #15, #T, RO ; RO = result = reserved operand 
O3DA 674 5 -0.0. RO will be copied t 
O3DA 675 3; signal mechanism vector (CHE SL _MCH_RO/R1) 
O3DA 676 3; so can be fixed up by any error 
O3DA 677 ; handler 
00000000'GF 01 FB O3DA 678 CALLS #1, G*MTHSS$SIGNAL 3 signal condition 
04 O3E1 679 RET ; return - RO restored from CHFSL_MCH_RO/R1 
O3E2 680 
O3E2 681 - END 


OTSSPOWDD 
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o 
ooo 
2Zz2 


Ai_ TABLE 0 
A2-TABLE 8 
eee = 
BAS = 
C 0 
C1 
C2 
C 
DEF INED 
AL_LOG2 
EXCEPTION_1 0 
EXCEPTION 2 0 
EXPLEN z 
EXPTAB 02 
EXP z= 
GET_BASE ¢ 
NDEX 0 
LOGLEN = 
TAB 02 
MTHSS$SIGNAL teereree x 00 
MTHSK_FLOOVEMAT teeeeeee xX 00 
MTHSK_FLOUNDMAT eteeeeee xX 00 
MTHSK_UNDEXP eeeerere §=6X 00 
S$POwDD 00000 RG 0 
OTSSPOWDR 00000218 RG 0 
OTSSPOWRD 00000232 RG 0 
VER 000003D2 R 0 
RETURN 00000 R 02 
SFSW_SAVE_PSW = 00000004 
TwOo_ASS 00000190 R 8s 
UNDEF INED siaiset tHe g 0 
f QOOOOSBF R 02 
Y_TIMES _W1_OVER 00000388 R 02 
+e ewe aaee aoe e neces ae + 
‘ ;_Psect synopsis ; 
PSECT name Allocation PSECT No. Attributes 
- ABS . 00000000 < 0.) 00 ¢ O.) NOPIC 
SABSS 00009000 = ( 0.) 0O1¢ 1.) NOPIC 
_OTSSCODE 000003E2 ( 994.) 02 ¢ 2.) PIC 
nasa annonce agai tacit Ae 
! ; Performance indicators ! 
Phase Page faults CPU Time Elapsed Time 
Initialization 4 0:00:00.1 0:00:00.63 
Command processing 1 :00:00.7 0: 8: 6. 
Pass 1 134 :00:02. :00: 8: 
Symbol table sort :00:00.05 :00:00.1 
ass 13 :00:01.51 :00:04.5 
Symbol table output 0.06 0:00.54 


11:28:08 EMTHRTL.SAcSorSPouDD.man;1 "29° 


AB LCL poe he NOEXE NORD NOWRT NOVEC BYTE 
ABS LCL NOSHR EXE RD WRT NOVEC BYTE 
RE LCL SHR EXE RD NOWRT NOVEC QUAD 
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OTS$ dD 
VAX-11 Macro Run Statistics -$ MTHRTL.SRCJOTSPOWDD.MAR; 1 
Psect synopsis output 


Crosecteterence output 8 B898:88:96 98:99 :90:98 


The york ing *st Limit was 1200 pages. 
12217 bytes (24 pages) of virtual memory were used to buffer the intermediate code. 
There were 10 pages of symbol table space allocated to hold 62 non-local and 2 local symbols. 


41 source Lines were read in Pass 1, producing 21 object records in Pass 2. 
9 pages of virtual memory were used to define 8 macros. 


oes 


! Macro Library statistics ! 


fe enw ewan eer eee ren oom oanes + 


Macro Library name Macros defined 
-$255S$DUA28: CSYSLIBISTARLET .MLB;2 4 


88 GETS were required to define 4 macros. 
There were no errors, warnings or information messages. 
MACRO/ENABLE=SUPPRESSION/D1I SABLE=(GLOBAL , TRACEBACK) /LIS=LIS$:OTSPOWDD/OBJ=0BJ$:OTSPOWDD MSRC$:MTHJACKET/UPDATE=(ENHS:MTHJACKET) +MSRC 
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